#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__1 = 1;


/* *********************************************************************** */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dlabcm_
   (integer *n, integer *nband, integer *nl, integer*nr, doublereal *a,
   doublereal *eigval, integer *lde, doublereal *eigvec, doublereal *atol,
   doublereal *artol, doublereal *bound, doublereal *atemp, doublereal *d__,
   doublereal *vtemp)
{
  /* System generated locals */
  integer a_dim1, a_offset, eigvec_dim1, eigvec_offset, i__1, i__2, i__3;
  doublereal d__1, d__2;

  /* Local variables */
  static logical flag__;
  static doublereal errb;
  /*extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
  integer *);*/
  static integer nval, numl;
  //extern doublereal dnrm2_(integer *, doublereal *, integer *);
  static integer i__, j, l, m;
  //extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
  // integer *);
  static doublereal sigma, resid;
  //extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *), daxpy_(integer *, doublereal *, 
  // doublereal *, integer *, doublereal *, integer *);
  static doublereal vnorm;
  //extern /* Subroutine */ int dlabfc_(integer *, integer *, doublereal *, 
  // doublereal *, integer *, integer *, doublereal *, integer *, 
  // integer *, doublereal *, doublereal *, doublereal *), dlabax_(
  // integer *, integer *, doublereal *, doublereal *, doublereal *), 
  // dlaran_(integer *, doublereal *);
  static doublereal rq;
  static integer numvec;
  static doublereal gap;


  /*  THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES */
  /*  FOR THE BNDEIG PACKAGE.  EIGENVALUES ARE COMPUTED BY */
  /*  A MODIFIED RAYLEIGH QUOTIENT ITERATION.  THE EIGENVALUE COUNT */
  /*  OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE */
  /*  THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO */
  /*  INSURE CONVERGENCE TO THE DESIRED EIGENVALUES. */

  /*  FORMAL PARAMETERS. */



  /*  LOCAL VARIABLES */



  /*  FUNCTIONS CALLED */


  /*  SUBROUTINES CALLED */

  /*     DLABAX, DLABFC, DLARAN, DAXPY, DCOPY, DSCAL */

  /*  REPLACE ZERO VECTORS BY RANDOM */

  /* Parameter adjustments */
  a_dim1 =  *nband;
  a_offset = 1+a_dim1 * 1;
  a -= a_offset;
  --eigval;
  eigvec_dim1 =  *lde;
  eigvec_offset = 1+eigvec_dim1 * 1;
  eigvec -= eigvec_offset;
  bound -= 3;
  --atemp;
  --d__;
  --vtemp;

  /* Function Body */
  nval =  *nr -  *nl + 1;
  flag__ = FALSE_;
  i__1 = nval;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (ddot_(n, &eigvec[i__ *eigvec_dim1 + 1], &c__1, &eigvec[i__ *eigvec_dim1
      + 1], &c__1) == 0.)
    {
      dlaran_(n, &eigvec[i__ *eigvec_dim1 + 1]);
    }
    /* L5: */
  }

  /*  LOOP OVER EIGENVALUES */

  sigma = bound[(nval + 1 << 1) + 2];
  i__1 = nval;
  for (j = 1; j <= i__1; ++j)
  {
    numl = j;

    /*  PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT */

    L10: dlabax_(n, nband, &a[a_offset], &eigvec[j *eigvec_dim1 + 1], &vtemp[1])
      ;
    vnorm = dnrm2_(n, &vtemp[1], &c__1);
    if (vnorm == 0.)
    {
      goto L20;
    }
    d__1 = 1. / vnorm;
    dscal_(n, &d__1, &vtemp[1], &c__1);
    d__1 = 1. / vnorm;
    dscal_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1);
    d__1 =  - sigma;
    daxpy_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1, &vtemp[1], &c__1);

    /*  LOOP OVER SHIFTS */

    /*  COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE */

    L20: vnorm = dnrm2_(n, &eigvec[j *eigvec_dim1 + 1], &c__1);
    if (vnorm != 0.)
    {
      goto L30;
    }
    dlaran_(n, &eigvec[j *eigvec_dim1 + 1]);
    goto L10;

    L30: rq = sigma + ddot_(n, &eigvec[j *eigvec_dim1 + 1], &c__1, &vtemp[1], 
      &c__1) / vnorm / vnorm;
    d__1 = sigma - rq;
    daxpy_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1, &vtemp[1], &c__1);
    /* Computing MAX */
    d__1 =  *atol, d__2 = dnrm2_(n, &vtemp[1], &c__1) / vnorm;
    resid = max(d__1, d__2);
    d__1 = (float)1. / vnorm;
    dscal_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1);

    /*  ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH */

    if (bound[(j + 1 << 1) + 2] - bound[(j + 1 << 1) + 1] <  *atol * 3.)
    {
      goto L300;
    }

    /*  COMPUTE MINIMAL ERROR BOUND */

    errb = resid;
    /* Computing MIN */
    d__1 = bound[(j + 2 << 1) + 1] - rq, d__2 = rq - bound[(j << 1) + 2];
    gap = min(d__1, d__2);
    if (gap > resid)
    {
      /* Computing MAX */
      d__1 =  *atol, d__2 = resid * resid / gap;
      errb = max(d__1, d__2);
    }

    /*  TENTATIVE NEW SHIFT */

    sigma = (bound[(j + 1 << 1) + 1] + bound[(j + 1 << 1) + 2]) *.5;

    /*  CHECK FOR TERMINALTION */

    if (resid >  *atol * 2.)
    {
      goto L40;
    }
    if (rq - errb > bound[(j << 1) + 2] && rq + errb < bound[(j + 2 << 1) + 1])
    {
      goto L310;
    }

    /*  RQ IS TO THE LEFT OF THE INTERVAL */

    L40: if (rq >= bound[(j + 1 << 1) + 1])
    {
      goto L50;
    }
    if (rq - errb > bound[(j << 1) + 2])
    {
      goto L100;
    }
    if (rq + errb < bound[(j + 1 << 1) + 1])
    {
      dlaran_(n, &eigvec[j *eigvec_dim1 + 1]);
    }
    goto L200;

    /*  RQ IS TO THE RIGHT OF THE INTERVAL */

    L50: if (rq <= bound[(j + 1 << 1) + 2])
    {
      goto L100;
    }
    if (rq + errb < bound[(j + 2 << 1) + 1])
    {
      goto L100;
    }

    /*  SAVE THE REJECTED VECTOR IF INDICATED */

    if (rq - errb <= bound[(j + 1 << 1) + 2])
    {
      goto L200;
    }
    i__2 = nval;
    for (i__ = j; i__ <= i__2; ++i__)
    {
      if (bound[(i__ + 1 << 1) + 2] > rq)
      {
        goto L70;
      }
      /* L60: */
    }
    goto L80;

    L70: dcopy_(n, &eigvec[j *eigvec_dim1 + 1], &c__1, &eigvec[i__ *eigvec_dim1
      + 1], &c__1);

    L80: dlaran_(n, &eigvec[j *eigvec_dim1 + 1]);
    goto L200;

    /*  PERTURB RQ TOWARD THE MIDDLE */

    L100: if (sigma < rq)
    {
      /* Computing MAX */
      d__1 = sigma, d__2 = rq - errb;
      sigma = max(d__1, d__2);
    }
    if (sigma > rq)
    {
      /* Computing MIN */
      d__1 = sigma, d__2 = rq + errb;
      sigma = min(d__1, d__2);
    }

    /*  FACTOR AND SOLVE */

    L200: i__2 = nval;
    for (i__ = j; i__ <= i__2; ++i__)
    {
      if (sigma < bound[(i__ + 1 << 1) + 1])
      {
        goto L220;
      }
      /* L210: */
    }
    i__ = nval + 1;
    L220: numvec = i__ - j;
    /* Computing MIN */
    i__2 = numvec, i__3 =  *nband + 2;
    numvec = min(i__2, i__3);
    if (resid <  *artol)
    {
      numvec = min(1, numvec);
    }
    dcopy_(n, &eigvec[j *eigvec_dim1 + 1], &c__1, &vtemp[1], &c__1);
    i__2 = (*nband << 1) - 1;
    dlabfc_(n, nband, &a[a_offset], &sigma, &numvec, lde, &eigvec[j
      *eigvec_dim1 + 1], &numl, &i__2, &atemp[1], &d__[1], atol);

    /*  PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW */

    if (numvec == 1)
    {
      goto L227;
    }
    l = numvec - 1;
    i__2 = l;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      m = j + i__;
      d__1 = 1. / vnorm;
      dscal_(n, &d__1, &eigvec[m *eigvec_dim1 + 1], &c__1);
      /* L225: */
    }

    /*  UPDATE INTERVALS */

    L227: numl = numl -  *nl + 1;
    if (numl >= 0)
    {
      bound[4] = min(bound[4], sigma);
    }
    i__2 = nval;
    for (i__ = j; i__ <= i__2; ++i__)
    {
      if (sigma < bound[(i__ + 1 << 1) + 1])
      {
        goto L20;
      }
      if (numl < i__)
      {
        bound[(i__ + 1 << 1) + 1] = sigma;
      }
      if (numl >= i__)
      {
        bound[(i__ + 1 << 1) + 2] = sigma;
      }
      /* L230: */
    }
    if (numl < nval + 1)
    {
      /* Computing MAX */
      d__1 = sigma, d__2 = bound[(nval + 2 << 1) + 1];
      bound[(nval + 2 << 1) + 1] = max(d__1, d__2);
    }
    goto L20;

    /*  ACCEPT AN EIGENPAIR */

    L300: dlaran_(n, &eigvec[j *eigvec_dim1 + 1]);
    flag__ = TRUE_;
    goto L310;

    L305: flag__ = FALSE_;
    rq = (bound[(j + 1 << 1) + 1] + bound[(j + 1 << 1) + 2]) *.5;
    i__2 = (*nband << 1) - 1;
    dlabfc_(n, nband, &a[a_offset], &rq, &numvec, lde, &eigvec[j *eigvec_dim1 +
      1], &numl, &i__2, &atemp[1], &d__[1], atol);
    vnorm = dnrm2_(n, &eigvec[j *eigvec_dim1 + 1], &c__1);
    if (vnorm != (float)0.)
    {
      d__1 = 1. / vnorm;
      dscal_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1);
    }

    /*  ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES */

    L310: eigval[j] = rq;
    if (j == 1)
    {
      goto L330;
    }
    m = j - 1;
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      d__1 =  - ddot_(n, &eigvec[i__ *eigvec_dim1 + 1], &c__1,
        &eigvec[j*eigvec_dim1 + 1], &c__1);
      daxpy_(n, &d__1, &eigvec[i__ *eigvec_dim1 + 1], &c__1,
        &eigvec[j*eigvec_dim1 + 1], &c__1);
      /* L320: */
    }
    L330: vnorm = dnrm2_(n, &eigvec[j *eigvec_dim1 + 1], &c__1);
    if (vnorm == 0.)
    {
      goto L305;
    }
    d__1 = 1. / vnorm;
    dscal_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1);

    /*   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE */

    if (flag__)
    {
      goto L305;
    }
    if (j == nval)
    {
      return 0;
    }
    m = j + 1;
    i__2 = nval;
    for (i__ = m; i__ <= i__2; ++i__)
    {
      d__1 =  - ddot_(n, &eigvec[j *eigvec_dim1 + 1], &c__1,
        &eigvec[i__*eigvec_dim1 + 1], &c__1);
      daxpy_(n, &d__1, &eigvec[j *eigvec_dim1 + 1], &c__1,
        &eigvec[i__*eigvec_dim1 + 1], &c__1);
      /* L340: */
    }
    /* L400: */
  }
  return 0;

  /* L500: */
  return 0;
} /* dlabcm_ */
